Dynamical Phase Transitions and Instabilities in Open Atomic Many-Body Systems 
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We discuss an open driven-dissipative many-body system, in which the competition of unitary 
Hamiltonian and dissipative Liouvillian dynamics leads to a nonequilibrium phase transition. It 
shares features of a quantum phase transition in that it is interaction driven, and of a classical 
phase transition, in that the ordered phase is continuously connected to a thermal state. Within a 
generalized Gutzwiller approach which includes the description of mixed state density matrices, we 
characterize the complete phase diagram and the critical behavior at the phase transition approached 
as a function of time. We find a novel fluctuation induced dynamical instability, which occurs at long 
wavelength as a consequence of a subtle dissipative renormalization effect on the speed of sound. 
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Experiments with cold atoms provide a unique setting 
to study nonequilibrium phenomena and dynamics, both 
in closed systems but also for (driven) open quantum 
dynamics. This relics on the ability to control the many- 
body dynamics and to prepare initial states far from the 
ground state. For closed systems we have seen a plethora 
of studies of quench dynamics [TJE], thermalization [3j|4], 
and transport [5], and also dynamical studies of crossing 
in a finite time quantum critical points in the spirit of the 
Kibble-Zurek mechanism [5J [TJ. On the other hand, sys- 
tems of cold atoms can be driven by external (light) fields 
and coupled to dissipative baths, thus realizing driven 
open quantum systems. As familiar e.g. from the quan- 
tum optics of the laser, the steady state of such a system 
(if it exists) is characterized by a dynamical equilibrium 
between pumping and dissipation, and can exhibit vari- 
ous nonequilibrium phases and phase transitions |5] as 
function of external control parameters. In the present 
work we will study such scenarios for quantum degener- 
ate gases. Our emphasis is on understanding quantum 
phases and dynamical phase transitions of cold atoms as 
an interacting many-body condensed matter system far 
from equilibrium. 

For a many-body system in thermodynamic equilib- 
rium the competition of two noncommuting parts of a 
microscopic Hamiltonian H = H\ + gH2 manifests it- 
self as a quantum phase transition (QPT), if the ground 
states for g <C g c and g 3> g c have different symme- 
tries |10j . For temperature T = the critical value g c 
then separates two distinct quantum phases, while for fi- 
nite temperature this defines a quantum critical region 
around g c in a T vs. g phase diagram. A seminal ex- 
ample in the context of cold atoms in optical lattices 
is the superfluid-Mott insulator transition in the Bosc- 
Hubbard (BH) model, with Hamiltonian 

ff = -J^&]&^-/i^n* + i[/^n^-l) , (1) 
with hi bosonic operators annihilating a particle on site 



I, fi£ = blbi number operators, J the hopping amplitude, 
and U the onsite interaction strength. For a given chem- 
ical potential p, chosen to fix a mean particle density n, 
the critical coupling strength g c = (U/ Jz) c separates a 
superfluid Jz 3> U from a Mott insulator regime Jz <C U 
(z the lattice coordination number). 

In contrast, we consider a nonequilibrium situation in 
which the competition of microscopic quantum mechani- 
cal operators results from an interplay of unitary (Hamil- 
tonian) and dissipative (Liouvillian) dynamics. We study 
a cold atom evolution described by a master equation for 
the many-body density operator 

d tP = -i[H,p]+C[p] , (2) 

£[p\ = 2 K ^ f 2c ii'P4i' ~ c u> c U'P - pc\ e cwj , 

where cw = (b\ + b\,)(bi — by) are Lindblad "jump 
operators" acting on adjacent sites The energy 

scale k is the dissipative rate. As shown in [TTj . such 
dissipative reservoir couplings are obtained in a setup 
where laser driven atoms are coupled to a phonon bath 
provided by a second condensate. For no interaction 
U = this dissipation drives the system to a dynami- 
cal equilibrium independent of the initial state [11] given 
by the pure many body state p ss = |BEC)(BEC| rep- 
resenting a Bose Einstein condensate. From an atomic 
physics point of view this is remarkable, as typical deco- 
herence mechanisms, such as spontaneous emission act- 
ing locally on lattice sites, will destroy long range or- 
der, whereas here the bath coupling is engineered to sup- 
press phase fluctuations. This can be easily understood 
in momentum space, where the annihilation part of cw 
reads ~ exp^A ))^ with A the reciprocal lattice 

directions and a the lattice constant, cw thus feature 
a (unique) dissipative zero mode at q = - a many- 
body "dark state" |BEC) ~ &q^ |vac) decoupled from 
the bath, into which the system is consequently driven 
for long wait times. The dynamics behind Eq. (I3J) can 



2 



thus be understood as a "dark state laser cooling" [T2] 
into a condensate, although in a many-body context. 

|BEC) is also an eigenstate of kinetic energy. In 
contrast, turning on an interaction measured by u — 
U I (4kz) provides a Hamiltonian term in (J3| which is 
incompatible with kinetic energy and dissipation. This 
competition leads to novel dynamical equilibria which 
cannot be understood as thermodynamic equilibrium 
states found from minimizing a free energy. They are 
summarized in the steady state phase diagram in Fig. 
[T] Most prominently, it features a strong coupling phase 
transition as a function of u. A first hallmark of the 
nonequilibrium nature of the system is this: The tran- 
sition shares features of a QPT in that it is interaction 
driven, and of a classical phase transition in that the or- 
dered phase terminates in a mixed state. This contrasts 
e.g. the well-known dissipation induced phase transition 
to a superconductor in Josephson junction arrays [15] . 
in which detailed balance guarantees that the system's 
state remains pure despite the suppression of phase fluc- 
tuations via the coupling to a zero temperature bath. 

Furthermore, we show the existence of a novel dynam- 
ical instability that covers an extensive domain of the 
phase diagram. Again, this is a nonequilibrium effect, 
since in equilibrium, finite momentum excitations carry 
positive kinetic energy ruling out dynamical instabilities. 
It persists at arbitrarily weak interaction parameters U n 
due to its fluctuation induced nature elucidated below. 
This is in marked contrast to the "classical" dynamical 
instabilities of condensates in boosted lattices [HI [15] or 
in exciton-polariton systems [16] . which are induced by 
external tuning of parameters beyond finite critical val- 
ues. 

Nonlinear mean field master equation. — To solve the 
master equation we developed a generalized Gutzwiller 
approach, expected to hold in sufficiently high spatial 
dimension, which allows to include density matrices cor- 
responding to mixed states. This is implemented by a 
product ansatz p = pt, with the reduced local density 
operators pe — Tr^g p. The equation of motion (EoM) 
reads 

d t pe = -i[hi, pi] + Ci[pc] , (3) 

with the local Hamiltonian hi = —JJ2(e>\e)((be')b\ + 
(b\,)bi) — png + \l}hg(ng — 1) reproducing the stan- 
dard form of the Gutzwiller mean field approxi- 
mation and a Liouvillian of the form £i[pi] = 
«E(fW Et : s=i^l2Alp t A?-A?Alp t -p e A?Ar}. The 
Liouvillian is constructed with the vector of operators 
A? = (1, b\, bi, fit) and the matrix of correlation functions 

T r £ ' s = o- r o- s TreA?-^A?- r) pt, for a = (-1,-1,1,1). 
The p-dependent correlation matrix makes the master 
equation nonlinear in pp. 

Dynamical quantum phase transition. — At U = a 
steady state solution of Eq. ^ is given by the pure 
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FIG. 1: (color online) Nonequilibrium phase diagram for the 
model in Eq. |3|. The solid lines indicate the border of the 
dynamical quantum phase transition from a condensed to a 
homogeneous thermal steady state. The dashed lines delimit 
the region where the condensed state is stable with respect to 
spatial fluctuations. The black (blue) lines are the numerical 
results corresponding to average density n = 1.0 (n = 0.1). 
The red line corresponds to the analytical results for n — 0.1. 

state p^ — | V E')( V I'| for any I together with the choice 
p = —Jz, where \*f/) is a coherent state of parameter 
ne for any phase 9 [IT]- In order to understand the 
effect of a finite interaction U, we apply the rotating- 
frame transformation V(U) = exp[i[/n^(n£ — l)i] to 
Eq. ([§. This removes the interaction term from the 
unitary evolution, but the annihilation operators become 
VbgV -1 = J2 m exp(imUt)\m}g(m\bg. The effect of a fi- 
nite U is thus to rotate the phase of each Fock state 
differently, leading to dephasing of the coherent state 

(c) 

p\ . Hence, for strong enough [/, off-diagonal order is 
suppressed completely and the density matrix becomes 
diagonal. In this case Eq. ^ reduces precisely to the 
master equation for a system of bosons coupled to a ther- 
mal reservoir with occupation n [17) . whose solution is 
a mixed diagonal thermal state pw. Interestingly, this 
state is thermal-like; however the role of the thermal bath 
is played by the system itself, being provided by the mean 
occupation of neighbouring sites. 

We substantiate the discussion above with the numeri- 
cal integration of the EoM ^ for a homogeneous system 
(we drop the index £). The system is initially in the co- 
herent state and the condensate fraction \ip\ 2 /n, where 
ijj = (£>), decreases in time depending on the value of the 
interaction strength U. The result is a continuous tran- 
sition from the coherent state to the thermal state 
p^\ shown in Fig. [2] for some typical parameters. The 
boundary between the thermal and the condensed phase 
with varying J, n is shown in Fig. [T] with solid lines. 

The transition is a smooth crossover for any finite time, 
but for t — > 00 a sharp nonanalytic point indicating a 
second order phase transition develops. In the universal 
vicinity of the critical point, 1/nt may be viewed as an 
irrelevant coupling in the sense of the renormalization 
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FIG. 2: Stroboscopic plot of the time evolution of the conden- 
sate fraction as a function of the interaction strength U, for 
J = 1.5 k and n = 1. For large times it converges to the lower 
thick solid line. The critical point is U c — 4.5 kz. Inset: Near 
critical evolution reflected by the logarithmic derivative of the 
order parameter ip(t), for J = 0, n = 1, and U < U c . The 
early exponential decay (x) is followed by a scaling regime 
(o) with exponent a ~ 0.5. The final exponential runaway 
(+) is due to a small deviation from the critical point. 

group. We may use this attractive irrelevant direction to 
extract the critical exponent a for the order parameter 
from the scaling solution \ip(t)\ oc (nt)~ a . In the inset 
of Fig. [2] we plot a(t) = dlog(ip)/d\og(l/t) and read 
off the critical exponent a = 0.5 in the scaling regime, 
which is an expected result given the mean field nature 
of the Gutzwiller ansatz. We emphasize that following 
the relaxation dynamics of the condensate fraction for 
critical system parameters gives an experimental handle 
for the measurement of a. 

Low-density limit. — In the low density limit n -C 1 
we obtain an analytical understanding of the time evo- 
lution based on the observation that the six correlation 
functions ip, (6|), (b\b 2 ), and complex conjugates, form a 
closed (nonlinear) subset which decouples from the a pri- 
ori infinite hierarchy of normal ordered correlation func- 
tions (b\ n bf). We first use this result to obtain ana- 
lytically the critical exponent a discussed above. For a 
homogeneous system with J = the EoMs read 

d t ip = i^ + {-iU + An)(tfb 2 ) -4ki/>*(& 2 ) , 
d t (tfb 2 ) = %nnilj + {-iU + in-&K){tfb 2 ) , 
d t (b 2 ) = (-iU + 2ifi-8n)(b 2 ) +8nip 2 . (4) 

The structure of the equations suggest that (b 2 ) decays 
much faster than the other correlations for U = U c , so 
that we may take dt(b 2 ) = and hence (b 2 ) oc ip 2 . At the 
critical point the two linear contributions to dtip vanish 
due to the zero mass eigenvalue at criticality and dtip oc 
Kip 2 ip* . It follows that ~ l/(4y / Kt) in agreement with 
the numerical result in Fig. [2] 

To study the interaction induced depletion of the con- 
densate fraction, it is convenient to use "connected" cor- 



relation functions, built with the fluctuation operator 
5b = b — -00- Here ipo is the constant value of the or- 
der parameter in the steady state, and (5b) = 0. From 
Q we obtain a closed linear system of EoMs, if ip is 
considered as a parameter, determined self-consistently 
from the identity n — (StfSb) + IV-'ol 2 - The value of the 
chemical potential is fixed to remove the driving terms 
in the equations for (5b), leading to /i = nU. This is an 
equilibrium condition similar to the vanishing of the mass 
of the Goldstone mode in a thermodynamic equilibrium 
system with spontaneous symmetry breaking. The solu- 
tion of the equations in steady state yields the condensate 
fraction 

l^ol 2 = 1 2u 2 (l + (3+u) 2 ) 

n 1 + u 2 + j(8u + 67 (1 + 2u 2 ) + 24fu + 8j 3 ) ' 

(5) 

with dimensionless variable j = J/ (4k). Eq. ([5| re- 
duces to the simple quadratic expression 1 — 2u 2 in the 
limit of zero hopping, with the critical point U C (J = 
0) = 4/cz/V2- The phase boundary, obtained by set- 
ting ip = in Eq. reads u c = j + y?l/2 + 2j 2 . 
Fig. [l] shows that these compact analytical results (solid 
red line) match the full numerics for small densities (solid 
blue line) , and also explain the qualitative features of the 
phase boundary for large densities. We note the absence 
of distinct commensurability effects for e.g. n — 1, tied to 
the fact that the interaction also plays the role of heating. 

Dynamical instability. — Numerically integrating the 
full EoM ([3]) with site-dependence (in one dimension for 
simplicity), we observe a dynamical instability, mani- 
festing itself at late times in a long wavelength density 
wave with growing amplitude. Numerical linearization of 
Eq. ([3]) around the homogeneous steady state allows to 
draw a phase border for the unstable phase (see Fig. ^J. 
The instability is cured by the increase of hopping J, 
which is associated to an operator compatible with dissi- 
pation k. Furthermore, we note that the thermal state is 
always dynamically stable against long wavelength per- 
turbations. 

The origin of this instability is intriguing and we 
discuss it analytically within the low-density limit 
introduced above. We linearize in time the EoM 
([3]), writing the generic connected correlation func- 
tion as (6 e ){t) = (6e)o + S(6 e )(t), where (6 e ) 
is evaluated on the homogeneous steady state of 
the system. The EoM for the time and space de- 
pendent fluctuations is then Fourier transformed, 
resulting in a 7x7 matrix evolution equation 
dtS^q = A/<5<I>q for the correlation functions <j> q = 
((Sb) v (Sb%, (StfSb)^ 2 )^ 2 )^ (5bW) v (6bWSb)J. 
We note that the fluctuation 5(5b)q (<5(<56 r ) q ) coincides 
with the fluctuation of the order parameter Sipq (5ip*). 
The full matrix M can be easily diagonalized numerically 
revealing the spectrum in Fig. [3] (we display only the 
relevant real part 7 corresponding to damping). The 
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FIG. 3: Real (dissipative) part of the spectrum 7 q from the 
analytical low density limit for J = 0, n = 0.1, and U = 1.0 k. 
The inset magnifies the parameter region with unstable modes 
(red solid line). The black solid line is the bare dissipative 
spectrum ft q . 

lowest-lying branch gives 7 q < in a small interval 
around the origin q — 0. This means that the correlation 
functions grow exponentially oc e 7 * in a range of low 
momenta, resulting e.g. in a long wavelength density 
wave. 

Due to the scale separation for q — > in the matrix 
M apparent from Fig. |3j we can apply second order per- 
turbation theory twice in a row to integrate out the fast 
modes 7 oc k and oc nn. We then obtain an effective low 
energy EoM for the fluctuations of the order parameter 
(<^/> q , £^>* q ), governed by a 2 x 2 matrix 

_ / Un + tq - iKq Un + QwiKq \ , . 

cff I — Un — 9urtK q — Un — e q — ift q J ' 

where e q = Jq 2 represents the kinetic contribution and 
Kq = 2{2n + l)«q 2 is the bare dissipative spectrum. The 
form of the EoM reflects the structure of the spatial fluc- 
tuations which are included in our approach, that may be 
understood as scattering off the mean fields in opposite 
directions. We note that a naive a priori restriction to 
the 2 x 2 set corresponding to the subset (Sipe, 5ipg) would 
be inconsistent, for example destroying the dark state 
property present in the correct solution M e gf. On the 
other hand, factorizing the correlation functions in the 
Liouvillian Cg yields a dissipative Gross-Pitaevski equa- 
tion but its linearization in time produces a matrix M c g 
without the fluctuation induced term ~ u and fails to 
describe the dynamical instability. Thus, in order to cor- 
rectly capture the physics of the instability at long wave- 
length q — > 0, the onsite quantum correlations renormal- 
izing M e g have to be properly taken into account. 

We can make the nature of the instability even more 
transparent from calculating the lowest eigenvalue of 



M c g, 7 q ~ ic|q| + Kq, with speed of sound c = 
yj2Un[J — 9t/n/(2z)]. If the hopping amplitude is 
smaller than the critical value J c = 9Un/(2z) the speed 
of sound turns imaginary and contributes to the dissi- 
pative real part of 7 q . The nonanalytic renormalization 
contribution ~ |q| always dominates the bare quadratic 
piece for low momenta, explaining the shape in the inset 
of Fig. [3] and rendering the system unstable. The linear 
slope of the stability border for small J and U is clearly 
visible from the numerical results in Fig. [T] In summary, 
the origin of the instability is traced back to a subtle 
renormalization effect of the speed of sound at low ener- 
gies, which in turn is due to an interplay of short time 
quantum and long wavelength classical fluctuations. 

Conclusion. — We have discussed the steady state phase 
diagram resulting from a competition of unitary Bose- 
Hubbard and dissipative dynamics with dark state. The 
features found in the present model are expected to be 
generic and representative for a whole class of noncqui- 
librium models discussed recently in the context of reser- 
voir engineering and dissipative preparation of given long 
range ordered entangled states of qubits or spins on a 
lattice [HJ \T§\ and paired fermions [TTJ [30] . In particu- 
lar, we emphasize the importance of a compatible energy 
term for the achievement of stability of driven-dissipative 
many-body systems in future experiments. 

We thank M. Hayn, A. Pelster, S. Kehrein, M. Mockel, 
and J. V. Porto for interesting discussions. This work was 
supported by the Austrian Science Foundation through 
SFB FOQUS, SCALA and by EU Networks. 



[1] P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 
(2006); C. Kollath, A. M. LSuchli, and E. Altman, ibid. 
98, 180601 (2007); A. Silva, ibid. 101, 120603 (2008); M. 
Mockel and S. Kehrein, ibid. 100, 175702 (2008). 

[2] M. Greiner, O. Mandel, T.W. Hansen, and I. Bloch, Na- 
ture 419, 51 (2002); B. Paredes et al, ibid. 429, 277 
(2004); L.E. Sadler et al, ibid. 443, 312 (2006). 

[3] M. Cramer, CM. Dawson, J. Eisert, and T.J. Osborne, 
Phys. Rev. Lett. 100, 030602 (2008); M. Rigol, V. Dun- 
jko, and M. Olshanii, Nature 452, 854 (2008); G. Roux, 
Phys. R ev. A 79 , 021608(R) (2009); L.C. Venuti and P. 
Zanardi, |arXiv:0912.3357l 

[4] T. Kinoshita, T. Wenger, and D.S. Weiss, Nature 440, 
900 (2006); S. Hofferberth et al, Nature Phys. 4, 489 
(2008). 

[5] S. Montangero, R. Fazio, P. Zoller, and G. Pupillo, Phys. 
Rev. A 79, 041602(R) (2009); J. Schachenmayer, G. 
Pupillo, and A.J. Daley, New J. Phys. 12, 025014 (2010). 

[6] K. Sengupta, S. Powell, and S. Sachdev, Phys. Rev. A 
69, 053616 (2004); W.H. Zurek, U. Dorner, and P. Zoller, 
Phys. Rev. Lett. 95, 105701 (2005); T. Prosen and I. Pi- 
zorn, ibid. 101, 105701 (2008); C. De G randi, V. Gritsev 
and A. Polkovn ikov, |arXiv:0909.5181| R.A. Barankov, 
larXiv:0910. 02551 

[7] C.N. Weiler et al, Nature 455, 948 (2008). 



5 



[8] S.A. Moskalenko and D.W. Snoke, Bose-Einstein Con- 
densation of Excitons and Biexcitons, Cambridge Univ. 
Press, Cambridge (2000); J. Keeling, F.M. Marchetti, 
M.H. Szymanska, and P.B. Littlewood, Semicond. Sci. 
Technol. 22, Rl (2007). 
[9] E.G. Dalla Torre, E. D ernier, T. Giamarchi, and E. Alt- 
man, [arXw0908. 0868 

[10] S. Sachdev, Quantum Phase Transitions, Cambridge 
Univ. Press, Cambridge (1999). 

[11] S. Diehl et al, Nature Phys. 4, 1073 (2008); B. Kraus et 
ai, Phys. Rev. A 78, 042307 (2008). 

[12] A. Aspect et al, Phys. Rev. Lett. 61, 826 (1988); M. 
Kasevich and S. Chu, ibid. 69, 1741 (1992). 

[13] A. Schmid, Phys. Rev. Lett. 51, 1506 (1983); S. 
Chakravarty, G.-L. Ingold, S. Kivelson, and A. Luther, 
ibid. 56, 2303 (1986); A. Kampf and G. Schon, Phys. Rev. 
B 36, 3651 (1987); S. Chakravarty, S. Kivelson, G.T. Zi- 
manyi, and B.I. Halperin, ibid. 35, 7256 (1987); R. Fazio 
and H. van der Zant, Phys. Rep. 355, 235 (2001). 



[14] B. Wu and Q. Niu, Phys. Rev. A 64, 061603(R) (2001); 

A. Smerzi, A. Trombettoni, P.G. Kevrekidis, and A.R. 

Bishop, Phys. Rev. Lett. 89, 170402 (2002); E. Altman 

et al, ibid. 95, 020402 (2005); A. Polkovnikov et al, 

Phys. Rev. A 71, 063613 (2005). 
[15] S. Burger et al, Phys. Rev. Lett. 86, 4447 (2001); M. 

Cristiani et al, Optics Express 12, 4 (2004); J. Mun et 

al, Phys. Rev. Lett. 99, 150604 (2007). 
[16] J. Kasprzak et al, Nature 443, 409 (2006); M. Wouters 

and I. Carusotto, | |arXiv:1001. 06601 
[17] C.W. Gardiner and P. Zoller, Quantum Noise , Springer- 

Verlag, Berlin (1999). 
[18] F. Verstraete, M.M. Wolf, and J.I. Cirac, Nature Phys. 

5, 633 (2009). 

[19] H. Weimer et al, Nature Phys., doi:10.1038/nphysl614 
(2010). 

[20] S. Diehl, W. Yi, A. J. Daley, P. Zoller, |arXrv:1007.3420| 
(2010). 



